Scenario

Within this chapter developmental plasticity was explored within Acanthochromis polyacanthus that were collected from two different regions (i.e., low-latitude, Cairns, and high-latitude, Mackay). Fish were held in common garden experiments at 28.5 C. Clutch data was collected along with parental morphometeric data to determine how fish from each population performed at 28.5 C, a temperature that was shown to produce similar absoluate aerobic scope performances in a previous study [link]. Once hatched offspring were placed into three different temperature treatments, 28.5, 30, and 31.5 C. After approximately 5-6 months offspring length and weight was measured, as well as CTmax and respiration during CTmax trials.

Load packages

library(tidyverse) # data manipulation
library(ggpubr) # figure arrangement 
library(brms) # Bayesian models
library(StanHeaders)# needed to run Bayesian models
library(rstan) # needed to run Bayesian models
library(standist) # needs to be installed 
library(bayesplot) # needed for MCMC diagnostics 
library(DHARMa) # model validation 
library(ggdist) # partial plots 
library(tidybayes) # partial plots 
library(broom.mixed) # model investigation
library(emmeans) # pairwise comparisons
library(rstanarm) # pairwise comparisons - need for emmeans 
library(ggpubr) # arranging figures 
library(mgcv) # GAM 
library(modelr) # plotting
library(segmented) # piece wise regression 
library(strucchange) #piece wise regression

Set working directory

knitr::opts_knit$set(root.dir=working.dir)

Import data

This data set has no point outliers however, in the next chunk 3 samples will be discarded due to poor data quailty

lat_resp_dat <- read_delim("C:/Users/jc527762/OneDrive - James Cook University/PhD dissertation/Data/Chapter4_Latitude/import_files/lat_resp_dat_no_outliers.txt", 
    delim = "\t", escape_double = FALSE, 
    col_types = cols(area = col_skip(), rate.a.spec = col_skip()), 
    trim_ws = TRUE)

Data manipulation

rows_of_data <- count(lat_resp_dat, sampleID)
lat_resp_dat2 <- lat_resp_dat |> 
  mutate(dev.temp = as.factor(dev.temp), 
         replicate = str_sub(sampleID, -1,-1), 
         population = factor(population)) |>
                                           # number of observations = 5758
  filter(sampleID != "56.CARL.137.28,5.1", # 5777 - 76 = 5682
         sampleID != "56.CARL.137.28,5.2", # 5701 - 64 = 5618
         sampleID != "60.LCKM.152.30.1"  # 5637 - 76 = 5542
  ) |> 
  filter(time_lag_sec >2001) |> # remove samples from first 5 cycles (i.e., first 33 minutes) 
  group_by(sampleID) |> 
  mutate(max_value_index = which.max(rate.output), 
         row_number = row_number(), 
         max.rate = max(rate.output), 
         low.rate = mean(rate.output[order(rate.output)[1:10]]), 
         net.rate = max.rate - low.rate) |> 
  filter(row_number <= max_value_index) |>
  ungroup() |> 
  mutate(region = factor(region), 
         dev.temp =factor(dev.temp), 
         level = as.factor(case_when(tank >= 1 & tank <= 199 ~ 1,
                           tank >= 200 & tank <= 299 ~ 2,
                           tank >= 300 & tank <= 399 ~ 3,
                           TRUE ~ NA_real_)), 
         female = factor(female)) |> 
  unite("dev.temp.region",c(dev.temp,region), remove=FALSE) |> 
  mutate(dev.temp.region = factor(dev.temp.region))

Exploratory data analysis

2nd order polynomial

eda1 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output)) + 
  geom_point(alpha=0.5, color = "black") + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 2, raw=TRUE), color = "red") + 
  theme_classic() +
  ggtitle("All data points - 2nd order polynomial")

eda2 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) + 
  geom_point(alpha=0.5) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 2, raw=TRUE), color = "red") + 
  facet_wrap(~ region) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Region - 2nd order polynomial")

eda3 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) + 
  geom_point(alpha=0.1) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 2, raw=TRUE)) + 
  facet_wrap(~ region) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Region - 2nd order polynomial") 

eda4 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) + 
  geom_point(alpha=0.1) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 2, raw=TRUE)) + 
  facet_wrap(~ dev.temp) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Region - 2nd order polynomial")

eda5 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=population)) + 
  geom_point(alpha=0.5) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 2, raw=TRUE), color = "red") + 
  facet_wrap(~ population) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Population - 2nd order polynomial")

eda6 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) + 
  geom_point(alpha=0.1) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 2, raw=TRUE)) + 
  facet_wrap(~ population) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Population - 2nd order polynomial")

eda.fig <- ggarrange(eda1,eda2,eda3,eda4,eda5,eda6,
          nrow = 6, 
          ncol=1); eda.fig

3rd order polynomial

eda1 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output)) + 
  geom_point(alpha=0.5, color = "black") + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 3, raw=TRUE), color = "red") + 
  theme_classic() +
  ggtitle("All data points - 2nd order polynomial")

eda2 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) + 
  geom_point(alpha=0.5) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 3, raw=TRUE), color = "red") + 
  facet_wrap(~ region) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Region - 2nd order polynomial")

eda3 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) + 
  geom_point(alpha=0.1) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 3, raw=TRUE)) + 
  facet_wrap(~ region) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Region - 2nd order polynomial") 

eda4 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) + 
  geom_point(alpha=0.1) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 3, raw=TRUE)) + 
  facet_wrap(~ dev.temp) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Region - 2nd order polynomial")

eda5 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=population)) + 
  geom_point(alpha=0.5) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 3, raw=TRUE), color = "red") + 
  facet_wrap(~ population) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Population - 2nd order polynomial")

eda6 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) + 
  geom_point(alpha=0.1) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 3, raw=TRUE)) + 
  facet_wrap(~ population) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Population - 2nd order polynomial")

eda.fig <- ggarrange(eda1,eda2,eda3,eda4,eda5,eda6,
          nrow = 6, 
          ncol=1); eda.fig

Fit model [random factors]

First we will place random factors only within out model and see if they do a good job explaining the variance within our model. We will also be include priors which will be based off of out length data.

Hypothesis test will include:

Null model
EGG_COUNT ~ 1 + (1| FEMALE) + (1| TANK)
EGG_COUNT ~ 1 + (1| FEMALE) + (1| TANK) + (1| LEVEL)
EGG_COUNT ~ 1 + (1| FEMALE) + (1| TANK) + (1| LEVEL) + (1| POPULATION)
EGG_COUNT ~ 1 + (1| FEMALE) + (1| TANK) + (1| LEVEL) + (1| POPULATION)+ (1| CLUTCH_ORDER)

lat_resp_dat2 |> 
  group_by(region,dev.temp) |> 
  summarise(mean(na.omit(rate.output)), 
            sd(na.omit(rate.output)))

## `summarise()` has grouped output by 'region'. You can override using the
## `.groups` argument.

Priors

rate.priors <- prior(normal(430, 0.25), class="Intercept") +  
  prior(student_t(3, 0, 195), class = "sigma")

Models

Null

f.model.null <- bf(rate.output ~ 1, 
                   family = gaussian()) 

model.null <- brm(f.model.null,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE),  
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95))

## Compiling Stan program...

## Start sampling

saveRDS(model.null, file="./03.CTmax_resp_data_analsyis_files/model.null.RDS")

Model1

f.model1 <- bf(rate.output ~ 1 + (1| female) + (1| tank), 
                          family=gaussian())

model1 <- brm(f.model1,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE),  
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95))

## Compiling Stan program...

## Start sampling

saveRDS(model1, file="./03.CTmax_resp_data_analsyis_files/model1.RDS")

Model2

f.model2 <- bf(rate.output ~ 1 + (1| female) + (1| tank) + (1| level), 
                          family=gaussian())

model2 <- brm(f.model2,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE),  
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95))

## Compiling Stan program...

## Start sampling

## Warning: There were 13 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.

## Warning: Examine the pairs() plot to diagnose sampling problems

saveRDS(model2, file="./03.CTmax_resp_data_analsyis_files/model2.RDS")

Model3

f.model3 <- bf(rate.output ~ 1 + (1| female) + (1| tank) + (1| level) + (1| population), 
               family=gaussian())

model3 <- brm(f.model3,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE),  
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95))

## Compiling Stan program...

## Start sampling

## Warning: There were 5 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.

## Warning: Examine the pairs() plot to diagnose sampling problems

saveRDS(model3, file="./03.CTmax_resp_data_analsyis_files/model3.RDS")

Model4

f.model4 <- bf(rate.output ~ 1 + (1| female) + (1 | tank) + (1| level) + (1| population)+ (1| clutch_order), 
                        family=gaussian())

model4 <- brm(f.model4,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE),  
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95))

## Compiling Stan program...

## Start sampling

## Warning: There were 3 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.

## Warning: Examine the pairs() plot to diagnose sampling problems

saveRDS(model4, file="./03.CTmax_resp_data_analsyis_files/model4.RDS")

LOO

LOO(model.null,model1, model2,model3,model4)

## Output of model 'model.null':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -31469.1  78.1
## p_loo         2.5   0.4
## looic     62938.2 156.2
## ------
## MCSE of elpd_loo is 0.0.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.0]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model1':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -30662.6  99.6
## p_loo        57.6   3.1
## looic     61325.1 199.1
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.1]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model2':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -30662.7  99.5
## p_loo        57.6   3.0
## looic     61325.4 198.9
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model3':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -30662.7  99.6
## p_loo        58.1   3.2
## looic     61325.5 199.2
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model4':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -30662.5  99.4
## p_loo        57.5   3.0
## looic     61325.0 198.8
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Model comparisons:
##            elpd_diff se_diff
## model4        0.0       0.0 
## model1       -0.1       0.4 
## model2       -0.2       0.4 
## model3       -0.2       0.4 
## model.null -806.6      48.2

Fit model [fixed factors - polynomials]

Prior investigation

lat_resp_dat2 |> 
  group_by(region,dev.temp) |> 
  summarise(mean(na.omit(rate.output)), 
            sd(na.omit(rate.output)))

## `summarise()` has grouped output by 'region'. You can override using the
## `.groups` argument.

Priors

rate.priors <- prior(normal(430, 0.25), class="Intercept") + 
  prior(normal(0, 40), class="b")
  prior(student_t(3, 0, 195), class = "sigma")

Models

Model1.p1

f.model1.p1 <- bf(rate.output ~ 1 + 
                         dev.temp*region + Temperature +
                         scale(mass, center=TRUE, scale=TRUE) + 
                         (1| female) + (1| tank) , 
                         family=gaussian()) 

model1.p1 <- brm(f.model1.p1,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE), 
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95))

## Compiling Stan program...

## Start sampling

saveRDS(model1.p1, file="./03.CTmax_resp_data_analsyis_files/model1.p1.RDS")

DHARMa

#--- distribution check ---#
pp_check(model1.p1, type = 'dens_overlay_grouped', ndraws=150, group="region")

#--- DHARMa checks ---#
preds <- posterior_predict(model1.p1, ndraws=250, summary=FALSE)
model1.p1_resids <- createDHARMa(simulatedResponse = t(preds), 
                              observedResponse = lat_resp_dat2$rate.output, 
                              fittedPredictedResponse = apply(preds, 2, mean), 
                              integerResponse = 'student')
plot(model1.p1_resids) ; testDispersion(model1.p1_resids)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.98596, p-value = 0.704
## alternative hypothesis: two.sided

MCMC diagnostics

stan_trace(model1.p1$fit)

## 'pars' not specified. Showing first 10 parameters by default.

stan_ac(model1.p1$fit)

## 'pars' not specified. Showing first 10 parameters by default.

stan_rhat(model1.p1$fit)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_ess(model1.p1$fit)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_dens(model1.p1$fit, separate_chains = TRUE)

## 'pars' not specified. Showing first 10 parameters by default.

model1.p2

f.model1.p2 <- bf(rate.output ~ 1 + 
                         dev.temp*region + poly(Temperature, 2) +
                         scale(mass, center=TRUE, scale=TRUE) + 
                         (1| female) + (1| tank) , 
                         family=gaussian()) 

model1.p2 <- brm(f.model1.p2,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE), 
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95))

## Compiling Stan program...

## Start sampling

saveRDS(model1.p2, file="./03.CTmax_resp_data_analsyis_files/model1.p2.RDS")

DHARMa

#--- distribution check ---#
pp_check(model1.p2, type = 'dens_overlay_grouped', ndraws=150, group="region")

#--- DHARMa checks ---#
preds <- posterior_predict(model1.p2, ndraws=250, summary=FALSE)
model1.p2_resids <- createDHARMa(simulatedResponse = t(preds), 
                              observedResponse = lat_resp_dat2$rate.output, 
                              fittedPredictedResponse = apply(preds, 2, mean), 
                              integerResponse = 'student')
plot(model1.p2_resids) ; testDispersion(model1.p2_resids)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.98971, p-value = 0.688
## alternative hypothesis: two.sided

MCMC diagnostics

stan_trace(model1.p2$fit)

## 'pars' not specified. Showing first 10 parameters by default.

stan_ac(model1.p2$fit)

## 'pars' not specified. Showing first 10 parameters by default.

stan_rhat(model1.p2$fit)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_ess(model1.p2$fit)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_dens(model1.p2$fit, separate_chains = TRUE)

## 'pars' not specified. Showing first 10 parameters by default.

model1.p3

f.model1.p3 <- bf(rate.output ~ 1 + 
                         dev.temp*region + poly(Temperature, 3) +
                         scale(mass, center=TRUE, scale=TRUE) + 
                         (1| female) + (1| tank) , 
                         family=gaussian()) 

model1.p3 <- brm(f.model1.p3,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE), 
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95))

## Compiling Stan program...

## Start sampling

saveRDS(model1.p3, file="./03.CTmax_resp_data_analsyis_files/model1.p3.RDS")

DHARMa

#--- distribution check ---#
pp_check(model1.p3, type = 'dens_overlay_grouped', ndraws=150, group="region")

#--- DHARMa checks ---#
preds <- posterior_predict(model1.p3, ndraws=250, summary=FALSE)
model1.p3_resids <- createDHARMa(simulatedResponse = t(preds), 
                              observedResponse = lat_resp_dat2$rate.output, 
                              fittedPredictedResponse = apply(preds, 2, mean), 
                              integerResponse = 'student')
plot(model1.p3_resids) ; testDispersion(model1.p3_resids)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.98512, p-value = 0.616
## alternative hypothesis: two.sided

MCMC diagnostics

stan_trace(model1.p3$fit)

## 'pars' not specified. Showing first 10 parameters by default.

stan_ac(model1.p3$fit)

## 'pars' not specified. Showing first 10 parameters by default.

stan_rhat(model1.p3$fit)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_ess(model1.p3$fit)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_dens(model1.p3$fit, separate_chains = TRUE)

## 'pars' not specified. Showing first 10 parameters by default.

LOO

LOO(model1.p1, model1.p2,model1.p3)

## Output of model 'model1.p1':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -29946.0  94.1
## p_loo        60.3   3.1
## looic     59892.1 188.2
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.1]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model1.p2':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -30525.1 100.7
## p_loo        58.9   3.2
## looic     61050.1 201.3
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.2]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model1.p3':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -30516.1 100.1
## p_loo        58.7   3.1
## looic     61032.1 200.2
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.2]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Model comparisons:
##           elpd_diff se_diff
## model1.p1    0.0       0.0 
## model1.p3 -570.0      31.3 
## model1.p2 -579.0      31.5

Fit model [fixed factors - GAM]

GAM model k=4

f.model1.gamk4 <- bf(rate.output ~ 1 + 
                         t2(dev.temp,region,Temperature, k=4, bs="fs", by=dev.temp.region) + 
                         scale(mass, center=TRUE, scale=TRUE) + 
                         (1| female) + (1| tank), 
                         family=gaussian()) 

model1.gamk4 <- brm(f.model1.gamk4,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE), 
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95))

## Compiling Stan program...

## Start sampling

## Warning: There were 3 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.

## Warning: There were 1 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded

## Warning: Examine the pairs() plot to diagnose sampling problems

saveRDS(model1.gamk4, file="./03.CTmax_resp_data_analsyis_files/model1.gamk4.RDS")

DHARMa

#--- distribution check ---#
pp_check(model1.gamk4, type = 'dens_overlay_grouped', ndraws=150, group="dev.temp.region")

#--- DHARMa checks ---#
preds <- posterior_predict(model1.gamk4, ndraws=250, summary=FALSE)
model1.gamk4_resids <- createDHARMa(simulatedResponse = t(preds), 
                              observedResponse = lat_resp_dat2$rate.output, 
                              fittedPredictedResponse = apply(preds, 2, mean), 
                              integerResponse = 'student')
plot(model1.gamk4_resids) ; testDispersion(model1.gamk4_resids)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.97661, p-value = 0.472
## alternative hypothesis: two.sided

MCMC diagnostics

stan_trace(model1.gamk4$fit)

## 'pars' not specified. Showing first 10 parameters by default.

stan_ac(model1.gamk4$fit)

## 'pars' not specified. Showing first 10 parameters by default.

stan_rhat(model1.gamk4$fit)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_ess(model1.gamk4$fit)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_dens(model1.gamk4$fit, separate_chains = TRUE)

## 'pars' not specified. Showing first 10 parameters by default.

### {-}

GAM model k=3

f.model1.gamk3 <- bf(rate.output ~ 1 + 
                         t2(dev.temp,region,Temperature, k=3, bs="fs", by=dev.temp.region) + 
                         scale(mass, center=TRUE, scale=TRUE) + 
                         (1| female) + (1| tank), 
                         family=gaussian()) 

model1.gamk3 <- brm(f.model1.gamk3,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE), 
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95))

## Compiling Stan program...

## Start sampling

## Warning: There were 3 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.

## Warning: Examine the pairs() plot to diagnose sampling problems

saveRDS(model1.gamk3, file="./03.CTmax_resp_data_analsyis_files/model1.gamk3.RDS")

DHARMa

#--- distribution check ---#
pp_check(model1.gamk3, type = 'dens_overlay_grouped', ndraws=150, group="region")

#--- DHARMa checks ---#
preds <- posterior_predict(model1.gamk3, ndraws=250, summary=FALSE)
model1.gamk3_resids <- createDHARMa(simulatedResponse = t(preds), 
                              observedResponse = lat_resp_dat2$rate.output, 
                              fittedPredictedResponse = apply(preds, 2, mean), 
                              integerResponse = 'student')
plot(model1.gamk3_resids) ; testDispersion(model1.gamk3_resids)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.97942, p-value = 0.448
## alternative hypothesis: two.sided

MCMC diagnostics

stan_trace(model1.gamk3$fit)

## 'pars' not specified. Showing first 10 parameters by default.

stan_ac(model1.gamk3$fit)

## 'pars' not specified. Showing first 10 parameters by default.

stan_rhat(model1.gamk3$fit)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_ess(model1.gamk3$fit)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_dens(model1.gamk3$fit, separate_chains = TRUE)

## 'pars' not specified. Showing first 10 parameters by default.

### {-}

AIC(model1.gamk3, model1.gamk4)

## Warning in AIC.default(model1.gamk3, model1.gamk4): models are not all fitted
## to the same number of observations

Partial effects plots

conditional_effects

model1.gamk4 |> 
  conditional_effects(spaghetti = TRUE, ndraws=200)

Model investigation

summary

summary(model1.gamk4)

## Warning: There were 3 divergent transitions after warmup. Increasing
## adapt_delta above 0.95 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: rate.output ~ 1 + t2(dev.temp, region, Temperature, k = 4, bs = "fs", by = dev.temp.region) + scale(mass, center = TRUE, scale = TRUE) + (1 | female) + (1 | tank) 
##    Data: lat_resp_dat2 (Number of observations: 4776) 
##   Draws: 2 chains, each with iter = 5000; warmup = 500; thin = 5;
##          total post-warmup draws = 1800
## 
## Smooth Terms: 
##                                                             Estimate Est.Error
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_1)      420.23    192.27
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_2)      409.00    203.52
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_1)   807.91    370.10
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_2)   838.69    338.93
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_1)      672.50    274.67
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_2)      685.97    275.46
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_1)   910.31    389.52
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_2)   989.69    425.40
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_1)      423.24    186.68
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_2)      454.59    194.73
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_1)   592.65    257.89
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_2)   543.84    236.79
##                                                             l-95% CI u-95% CI
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_1)      204.95   935.83
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_2)      166.01   941.29
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_1)   392.88  1758.28
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_2)   413.66  1700.17
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_1)      350.27  1420.75
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_2)      338.90  1390.80
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_1)   457.39  1876.23
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_2)   482.79  2114.40
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_1)      202.85   861.71
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_2)      217.51   970.44
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_1)   291.71  1249.28
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_2)   264.18  1129.12
##                                                             Rhat Bulk_ESS
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_1)    1.00     1819
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_2)    1.00     1648
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_1) 1.00     1630
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_2) 1.00     1600
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_1)    1.00     1741
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_2)    1.00     1751
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_1) 1.00     1776
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_2) 1.00     1760
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_1)    1.00     1883
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_2)    1.00     1700
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_1) 1.00     1599
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_2) 1.00     1583
##                                                             Tail_ESS
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_1)        1829
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_2)        1746
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_1)     1529
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_2)     1662
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_1)        1761
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_2)        1706
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_1)     1714
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_2)     1744
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_1)        1911
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_2)        1708
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_1)     1606
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_2)     1563
## 
## Group-Level Effects: 
## ~female (Number of levels: 15) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)   118.65     31.65    65.99   189.82 1.00     1128     1388
## 
## ~tank (Number of levels: 52) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)    82.36     10.87    64.11   107.11 1.00     1276     1524
## 
## Population-Level Effects: 
##                                  Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept                          430.00      0.24   429.54   430.47 1.00
## scalemasscenterEQTRUEscaleEQTRUE    19.55      3.77    12.02    26.65 1.00
##                                  Bulk_ESS Tail_ESS
## Intercept                            1801     1707
## scalemasscenterEQTRUEscaleEQTRUE     1683     1613
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma   114.35      1.19   111.99   116.71 1.00     1771     1784
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

#car::Anova(model1.geo.gamma)

R2

model1.gamk4 |> bayes_R2(summary = FALSE) |> median_hdci()

tidyMCMC

tidyMCMC(model1.gamk4, estimate.method='median', conf.int=TRUE, conf.method='HPDinterval')

gather draws

#model1.re.wo |> get_variables()
model1.gamk4 |> gather_draws(`b_.*|sigma`, regex =TRUE) |> 
  median_hdci()

bayesplot

model1.gamk4 |> mcmc_plot(type='intervals')

emmeans - pariwise

model1.gamk4 |> emmeans(pairwise ~ region*dev.temp, type="response") |> pairs(by="dev.temp") |> summary(infer=TRUE)

## NOTE: A nesting structure was detected in the fitted model:
##     dev.temp.region %in% (dev.temp*region)

model1.gamk4 |> emmeans(pairwise ~ region*dev.temp, type="response") |> pairs(by="region") |> confint()

## NOTE: A nesting structure was detected in the fitted model:
##     dev.temp.region %in% (dev.temp*region)

model1.gamk4 |> emmeans(pairwise ~ region*dev.temp, type="response") |> regrid()

## NOTE: A nesting structure was detected in the fitted model:
##     dev.temp.region %in% (dev.temp*region)

##  region  dev.temp emmean lower.HPD upper.HPD
##  core    28          444       367       511
##  leading 28          419       303       511
##  core    30          467       399       538
##  leading 30          399       295       512
##  core    31          477       400       547
##  leading 31          485       393       577
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95

model1.gamk4 |> emmeans(pairwise ~ region*dev.temp, type="response") |> pairs(by="region") |> confint()

## NOTE: A nesting structure was detected in the fitted model:
##     dev.temp.region %in% (dev.temp*region)

emtrenns - pariwise

emtrends(model1.gamk4, specs=pairwise~dev.temp*region, var = "Temperature")

## NOTE: A nesting structure was detected in the fitted model:
##     dev.temp.region %in% (dev.temp*region)

## $emtrends
##  dev.temp region  Temperature.trend lower.HPD upper.HPD
##  28       core               44.872     31.78      57.1
##  30       core                2.862     -8.80      15.7
##  31       core               24.090     13.29      35.6
##  28       leading            -0.815    -14.74      12.8
##  30       leading             1.168    -13.99      16.4
##  31       leading             7.000     -5.06      19.1
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95 
## 
## $contrasts
##  contrast                                estimate lower.HPD upper.HPD
##  dev.temp28 core - dev.temp30 core          42.04    24.173     60.14
##  dev.temp28 core - dev.temp31 core          20.95     4.017     38.19
##  dev.temp28 core - dev.temp28 leading       45.62    26.541     63.58
##  dev.temp28 core - dev.temp30 leading       43.94    25.948     64.08
##  dev.temp28 core - dev.temp31 leading       38.00    20.584     55.40
##  dev.temp30 core - dev.temp31 core         -21.32   -37.216     -3.77
##  dev.temp30 core - dev.temp28 leading        3.65   -14.587     23.10
##  dev.temp30 core - dev.temp30 leading        1.71   -16.474     22.04
##  dev.temp30 core - dev.temp31 leading       -4.18   -20.285     13.35
##  dev.temp31 core - dev.temp28 leading       24.81     6.454     42.12
##  dev.temp31 core - dev.temp30 leading       23.23     5.496     42.63
##  dev.temp31 core - dev.temp31 leading       17.00    -0.205     32.99
##  dev.temp28 leading - dev.temp30 leading    -1.55   -20.579     19.27
##  dev.temp28 leading - dev.temp31 leading    -7.62   -25.674     10.96
##  dev.temp30 leading - dev.temp31 leading    -6.18   -25.762     12.80
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95

probabilities

mtsqst <- model1.gamk4 |> emmeans(pairwise ~ region*dev.temp)

## NOTE: A nesting structure was detected in the fitted model:
##     dev.temp.region %in% (dev.temp*region)

mtsqrt2 <- mtsqst$contrasts |> gather_emmeans_draws()
mtsqrt2 %>% group_by(contrast) %>% dplyr::summarise(Prob = sum(.value>0)/n())

summary figure

resp.plot <- lat_resp_dat2 |> 
  group_by(region, dev.temp,dev.temp.region) |> 
  data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature), 
                          to=max(lat_resp_dat2$Temperature), 
                          length.out=100), 
            mass=0.00148299) |> 
  add_fitted_draws(model1.gamk4, 
                   n=300, 
                   re_formula=NA) |>  
  unite("exp_group", c("region","dev.temp"), remove=FALSE) |> 
  ggplot(aes(x=Temperature, color=exp_group)) + 
  geom_line(aes(y=.value, group=paste(exp_group, .draw)), alpha=0.05) + 
  geom_smooth(aes(y=.value, color=exp_group), se=FALSE, size=1.5) + 
  theme_classic() + 
  facet_wrap(~dev.temp); resp.plot

resp.plot2 <- lat_resp_dat2 |> 
  group_by(region, dev.temp,dev.temp.region) |> 
  data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature), 
                          to=max(lat_resp_dat2$Temperature), 
                          length.out=100), 
            mass=0.00148299) |> 
  add_fitted_draws(model1.gamk4, 
                   n=300, 
                   re_formula=NA) |>  
  unite("exp_group", c("region","dev.temp"), remove=FALSE) |> 
  ggplot(aes(x=Temperature, color=exp_group)) + 
  geom_line(aes(y=.value, group=paste(exp_group, .draw)), alpha=0.05) + 
  geom_smooth(aes(y=.value, color=exp_group), se=FALSE, size=1.5) + 
  theme_classic() + 
  facet_wrap(~region); resp.plot2

resp.plot3 <- lat_resp_dat2 |> 
  group_by(region, dev.temp,dev.temp.region) |> 
  data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature), 
                          to=max(lat_resp_dat2$Temperature), 
                          length.out=100), 
            mass=0.00148299) |> 
  add_fitted_draws(model1.gamk4, 
                   n=300, 
                   re_formula=NA) |>  
  unite("exp_group", c("region","dev.temp"), remove=FALSE) |> 
  ggplot(aes(x=Temperature, color=exp_group)) + 
  geom_line(aes(y=.value, group=paste(exp_group, .draw)), alpha=0.05) + 
  geom_smooth(aes(y=.value, color=exp_group), se=FALSE, size=1.5) + 
  theme_classic(); resp.plot3

Piecewise linear regression

Low-latitude - 28.5C

posterior_draws <- lat_resp_dat2 |> 
  group_by(region, dev.temp,dev.temp.region) |> 
  data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature), 
                          to=max(lat_resp_dat2$Temperature), 
                          length.out=100), 
            mass=0.00148299) |> 
  add_linpred_draws(model1.gamk4, 
                   n=300, 
                   re_formula=NA) |>  
  unite("exp_group", c("region","dev.temp"), remove=FALSE) |> 
  filter(dev.temp.region == "28_core")  

spline_fit <- smooth.spline(posterior_draws$Temperature, posterior_draws$.linpred, nknots = 4) 

ggplot(posterior_draws, aes(x=Temperature, y=.linpred)) + 
  geom_point(color="black", alpha=0.05, size=1) + 
  geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y, 
                                        color="red"), 
              size=2)

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

strucchange analysis

spline.data <- as.data.frame(x= spline_fit$x) |> 
  mutate(y = spline_fit$y)|> 
  rename(x = `spline_fit$x`)

plot(spline.data)

fs.data <- Fstats(spline.data$y~1) 
bp.data <- breakpoints(spline.data$y~1)
summary(bp.data)

## 
##   Optimal (m+1)-segment partition: 
## 
## Call:
## breakpoints.formula(formula = spline.data$y ~ 1)
## 
## Breakpoints at observation number:
##                       
## m = 1            70   
## m = 2         60    81
## m = 3         54 70 85
## m = 4      40 55 70 85
## m = 5   15 40 55 70 85
## 
## Corresponding to breakdates:
##                               
## m = 1                 0.7     
## m = 2            0.6      0.81
## m = 3            0.54 0.7 0.85
## m = 4        0.4 0.55 0.7 0.85
## m = 5   0.15 0.4 0.55 0.7 0.85
## 
## Fit:
##                                                                
## m   0         1         2         3         4         5        
## RSS 1170431.8  216925.0   88460.8   49097.0   44088.4   41482.4
## BIC    1229.8    1070.4     989.9     940.3     938.7     941.8

plot(bp.data)

breakpoint.values=c(54,70,85); spline.data[breakpoint.values, c("x")]

## [1] 33.86963 35.04506 36.14703

segement analysis

Use out from strucchange summary to number of breakpoints as well as approximate location of breakpoints

my.lm <- lm(data=spline.data, y ~ x) 

spline.data[breakpoint.values, c("x")]

## [1] 33.86963 35.04506 36.14703

my.seg <- segmented(my.lm, 
                    seg.Z = ~x, 
                    psi=list(x=spline.data[breakpoint.values, c("x")]))
#summary(my.seg) 

my.seg$psi ; slope(my.seg)

##         Initial     Est.     St.Err
## psi1.x 33.86963 32.14501 0.02701804
## psi2.x 35.04506 33.52461 0.02942867
## psi3.x 36.14703 34.93830 0.03209560

## $x
##           Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 -15.086 0.53064 -28.429   -16.139   -14.032
## slope2  25.029 1.05370  23.754    22.937    27.122
## slope3  65.828 1.05370  62.474    63.736    67.921
## slope4  99.259 0.48165 206.080    98.303   100.220

my.fitted <- fitted(my.seg)
my.model <- data.frame(temperature = spline.data$x, rate.output= my.fitted)

ggplot(my.model, aes(x=temperature, y=rate.output)) + geom_line()

ggplot(posterior_draws, aes(x=Temperature, y=.linpred, color=.linpred)) +
  geom_point(alpha=0.05, size=1) + 
  geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y), color="red", 
              size=3) + ylab("rate.output")+
  geom_vline(xintercept = my.seg$psi[1,2]) + 
   geom_vline(xintercept = my.seg$psi[2,2])+ 
  geom_vline(xintercept = my.seg$psi[3,2])+ 
   geom_vline(xintercept = my.seg$psi[1,1], color="blue") + 
  geom_vline(xintercept = my.seg$psi[2,1], color="blue") + 
  geom_vline(xintercept = my.seg$psi[3,1], color="blue") + 
  geom_line(data=my.model, aes(x=temperature, y=rate.output), color="black", size=2) +
  theme_classic()

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

Low-latitude - 30C

posterior_draws <- lat_resp_dat2 |> 
  group_by(region, dev.temp,dev.temp.region) |> 
  data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature), 
                          to=max(lat_resp_dat2$Temperature), 
                          length.out=100), 
            mass=0.00148299) |> 
  add_linpred_draws(model1.gamk4, 
                   n=300, 
                   re_formula=NA) |>  
  unite("exp_group", c("region","dev.temp"), remove=FALSE) |> 
  filter(dev.temp.region == "30_core")  

spline_fit <- smooth.spline(posterior_draws$Temperature, posterior_draws$.linpred, nknots = 4) 

ggplot(posterior_draws, aes(x=Temperature, y=.linpred)) + 
  geom_point(color="black", alpha=0.05, size=1) + 
  geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y, 
                                        color="red"), 
              size=2)

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

strucchange analysis

spline.data <- as.data.frame(x= spline_fit$x) |> 
  mutate(y = spline_fit$y)|> 
  rename(x = `spline_fit$x`)
plot(spline.data)

fs.data <- Fstats(spline.data$y~1) 
bp.data <- breakpoints(spline.data$y~1)
summary(bp.data)

## 
##   Optimal (m+1)-segment partition: 
## 
## Call:
## breakpoints.formula(formula = spline.data$y ~ 1)
## 
## Breakpoints at observation number:
##                       
## m = 1               78
## m = 2            70 85
## m = 3         55 70 85
## m = 4   15    55 70 85
## m = 5   15 30 55 70 85
## 
## Corresponding to breakdates:
##                               
## m = 1                     0.78
## m = 2                 0.7 0.85
## m = 3            0.55 0.7 0.85
## m = 4   0.15     0.55 0.7 0.85
## m = 5   0.15 0.3 0.55 0.7 0.85
## 
## Fit:
##                                                                
## m   0         1         2         3         4         5        
## RSS 1289086.6  223458.9   95112.0   79870.4   72776.6   72689.6
## BIC    1239.4    1073.4     997.2     988.9     988.8     997.9

plot(bp.data)

breakpoint.values=c(54,70,85); spline.data[breakpoint.values, c("x")]

## [1] 33.86963 35.04506 36.14703

segement analysis

Use out from strucchange summary to number of breakpoints as well as approximate location of breakpoints

my.lm <- lm(data=spline.data, y ~ x) 

spline.data[breakpoint.values, c("x")]

## [1] 33.86963 35.04506 36.14703

my.seg <- segmented(my.lm, 
                    seg.Z = ~x, 
                    psi=list(x=spline.data[breakpoint.values, c("x")]))
#summary(my.seg) 

my.seg$psi ; slope(my.seg)

##         Initial     Est.     St.Err
## psi1.x 33.86963 31.47574 0.05887590
## psi2.x 35.04506 34.20755 0.02233563
## psi3.x 36.14703 35.37701 0.02265965

## $x
##             Est. St.Err.  t value CI(95%).l CI(95%).u
## slope1  26.77600 1.37500  19.4730  24.04500    29.507
## slope2   0.81219 0.58750   1.3825  -0.35463     1.979
## slope3  75.59800 2.06930  36.5340  71.48800    79.708
## slope4 154.35000 0.99772 154.7000 152.37000   156.330

my.fitted <- fitted(my.seg)
my.model <- data.frame(temperature = spline.data$x, rate.output= my.fitted)

ggplot(my.model, aes(x=temperature, y=rate.output)) + geom_line()

ggplot(posterior_draws, aes(x=Temperature, y=.linpred, color=.linpred)) +
  geom_point(alpha=0.05, size=1) + 
  geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y), color="red", 
              size=3) + ylab("rate.output")+
  geom_vline(xintercept = my.seg$psi[1,2]) + 
   geom_vline(xintercept = my.seg$psi[2,2])+ 
  geom_vline(xintercept = my.seg$psi[3,2])+ 
   geom_vline(xintercept = my.seg$psi[1,1], color="blue") + 
  geom_vline(xintercept = my.seg$psi[2,1], color="blue") + 
  geom_vline(xintercept = my.seg$psi[3,1], color="blue") + 
  geom_line(data=my.model, aes(x=temperature, y=rate.output), color="black", size=2) +
  theme_classic()

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

Low-latitude - 31C

posterior_draws <- lat_resp_dat2 |> 
  group_by(region, dev.temp,dev.temp.region) |> 
  data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature), 
                          to=max(lat_resp_dat2$Temperature), 
                          length.out=100), 
            mass=0.00148299) |> 
  add_linpred_draws(model1.gamk4, 
                   n=300, 
                   re_formula=NA) |>  
  unite("exp_group", c("region","dev.temp"), remove=FALSE) |> 
  filter(dev.temp.region == "31_core")  

spline_fit <- smooth.spline(posterior_draws$Temperature, posterior_draws$.linpred, nknots = 4) 

ggplot(posterior_draws, aes(x=Temperature, y=.linpred)) + 
  geom_point(color="black", alpha=0.05, size=1) + 
  geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y, 
                                        color="red"), 
              size=2)

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

strucchange analysis

spline.data <- as.data.frame(x= spline_fit$x) |> 
  mutate(y = spline_fit$y)|> 
  rename(x = `spline_fit$x`)
plot(spline.data)

fs.data <- Fstats(spline.data$y~1) 
bp.data <- breakpoints(spline.data$y~1)
summary(bp.data)

## 
##   Optimal (m+1)-segment partition: 
## 
## Call:
## breakpoints.formula(formula = spline.data$y ~ 1)
## 
## Breakpoints at observation number:
##                       
## m = 1            74   
## m = 2            64 83
## m = 3         54 70 85
## m = 4      40 55 70 85
## m = 5   22 40 55 70 85
## 
## Corresponding to breakdates:
##                                
## m = 1                 0.74     
## m = 2                 0.64 0.83
## m = 3            0.54 0.7  0.85
## m = 4        0.4 0.55 0.7  0.85
## m = 5   0.22 0.4 0.55 0.7  0.85
## 
## Fit:
##                                                          
## m   0        1        2        3        4        5       
## RSS 951373.7 176567.4  73635.2  44061.0  39736.9  39265.8
## BIC   1209.0   1049.8    971.6    929.4    928.3    936.3

plot(bp.data)

breakpoint.values=c(55,70,85); spline.data[breakpoint.values, c("x")]

## [1] 33.94309 35.04506 36.14703

segement analysis

Use out from strucchange summary to number of breakpoints as well as approximate location of breakpoints

my.lm <- lm(data=spline.data, y ~ x) 

spline.data[breakpoint.values, c("x")]

## [1] 33.94309 35.04506 36.14703

my.seg <- segmented(my.lm, 
                    seg.Z = ~x, 
                    psi=list(x=spline.data[breakpoint.values, c("x")]))
#summary(my.seg) 

my.seg$psi ; slope(my.seg)

##         Initial     Est.     St.Err
## psi1.x 33.94309 33.16127 0.02769747
## psi2.x 35.04506 34.42367 0.02225931
## psi3.x 36.14703 35.50138 0.02594152

## $x
##            Est. St.Err. t value CI(95%).l CI(95%).u
## slope1   5.4424 0.21502  25.311    5.0153    5.8694
## slope2  33.1200 0.89665  36.937   31.3390   34.9010
## slope3  75.5130 1.08240  69.767   73.3630   77.6630
## slope4 108.9600 0.53408 204.020  107.9000  110.0300

my.fitted <- fitted(my.seg)
my.model <- data.frame(temperature = spline.data$x, rate.output= my.fitted)

ggplot(my.model, aes(x=temperature, y=rate.output)) + geom_line()

ggplot(posterior_draws, aes(x=Temperature, y=.linpred, color=.linpred)) +
  geom_point(alpha=0.05, size=1) + 
  geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y), color="red", 
              size=3) + ylab("rate.output")+
  geom_vline(xintercept = my.seg$psi[1,2]) + 
   geom_vline(xintercept = my.seg$psi[2,2])+ 
  geom_vline(xintercept = my.seg$psi[3,2])+ 
   geom_vline(xintercept = my.seg$psi[1,1], color="blue") + 
  geom_vline(xintercept = my.seg$psi[2,1], color="blue") + 
  geom_vline(xintercept = my.seg$psi[3,1], color="blue") + 
  geom_line(data=my.model, aes(x=temperature, y=rate.output), color="black", size=2) +
  theme_classic()

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

High-latitude - 28C

posterior_draws <- lat_resp_dat2 |> 
  group_by(region, dev.temp,dev.temp.region) |> 
  data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature), 
                          to=max(lat_resp_dat2$Temperature), 
                          length.out=100), 
            mass=0.00148299) |> 
  add_linpred_draws(model1.gamk4, 
                   n=300, 
                   re_formula=NA) |>  
  unite("exp_group", c("region","dev.temp"), remove=FALSE) |> 
  filter(dev.temp.region == "28_leading")  

spline_fit <- smooth.spline(posterior_draws$Temperature, posterior_draws$.linpred, nknots = 4) 

ggplot(posterior_draws, aes(x=Temperature, y=.linpred)) + 
  geom_point(color="black", alpha=0.05, size=1) + 
  geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y, 
                                        color="red"), 
              size=2)

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

strucchange analysis

spline.data <- as.data.frame(x= spline_fit$x) |> 
  mutate(y = spline_fit$y)|> 
  rename(x = `spline_fit$x`)
plot(spline.data)

fs.data <- Fstats(spline.data$y~1) 
bp.data <- breakpoints(spline.data$y~1)
summary(bp.data)

## 
##   Optimal (m+1)-segment partition: 
## 
## Call:
## breakpoints.formula(formula = spline.data$y ~ 1)
## 
## Breakpoints at observation number:
##                       
## m = 1               79
## m = 2            70 85
## m = 3         55 70 85
## m = 4   15    55 70 85
## m = 5   15 38 55 70 85
## 
## Corresponding to breakdates:
##                                
## m = 1                      0.79
## m = 2                  0.7 0.85
## m = 3             0.55 0.7 0.85
## m = 4   0.15      0.55 0.7 0.85
## m = 5   0.15 0.38 0.55 0.7 0.85
## 
## Fit:
##                                                    
## m   0       1       2       3       4       5      
## RSS 1612756  264488  108644   99987   98860   97905
## BIC    1262    1090    1010    1011    1019    1028

plot(bp.data)

breakpoint.values=c(55,70,85); spline.data[breakpoint.values, c("x")]

## [1] 33.94309 35.04506 36.14703

segement analysis

Use out from strucchange summary to number of breakpoints as well as approximate location of breakpoints

my.lm <- lm(data=spline.data, y ~ x) 

spline.data[breakpoint.values, c("x")]

## [1] 33.94309 35.04506 36.14703

my.seg <- segmented(my.lm, 
                    seg.Z = ~x, 
                    psi=list(x=spline.data[breakpoint.values, c("x")]))
#summary(my.seg) 

my.seg$psi ; slope(my.seg)

##         Initial     Est.     St.Err
## psi1.x 33.94309 34.20914 0.05958793
## psi2.x 35.04506 34.93752 0.07430573
## psi3.x 36.14703 35.75777 0.08572172

## $x
##           Est. St.Err. t value CI(95%).l CI(95%).u
## slope1   1.386 0.61625  2.2490   0.16205    2.6099
## slope2  66.933 8.64990  7.7379  49.75300   84.1120
## slope3 135.080 7.49110 18.0320 120.20000  149.9600
## slope4 186.000 2.83140 65.6920 180.37000  191.6200

my.fitted <- fitted(my.seg)
my.model <- data.frame(temperature = spline.data$x, rate.output= my.fitted)

ggplot(my.model, aes(x=temperature, y=rate.output)) + geom_line()

ggplot(posterior_draws, aes(x=Temperature, y=.linpred, color=.linpred)) +
  geom_point(alpha=0.05, size=1) + 
  geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y), color="red", 
              size=3) + ylab("rate.output")+
  geom_vline(xintercept = my.seg$psi[1,2]) + 
   geom_vline(xintercept = my.seg$psi[2,2])+ 
  geom_vline(xintercept = my.seg$psi[3,2])+ 
   geom_vline(xintercept = my.seg$psi[1,1], color="blue") + 
  geom_vline(xintercept = my.seg$psi[2,1], color="blue") + 
  geom_vline(xintercept = my.seg$psi[3,1], color="blue") + 
  geom_line(data=my.model, aes(x=temperature, y=rate.output), color="black", size=2) +
  theme_classic()

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

High-latitude - 30C

posterior_draws <- lat_resp_dat2 |> 
  group_by(region, dev.temp,dev.temp.region) |> 
  data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature), 
                          to=max(lat_resp_dat2$Temperature), 
                          length.out=100), 
            mass=0.00148299) |> 
  add_linpred_draws(model1.gamk4, 
                   n=300, 
                   re_formula=NA) |>  
  unite("exp_group", c("region","dev.temp"), remove=FALSE) |> 
  filter(dev.temp.region == "30_leading")  

spline_fit <- smooth.spline(posterior_draws$Temperature, posterior_draws$.linpred, nknots = 4) 

ggplot(posterior_draws, aes(x=Temperature, y=.linpred)) + 
  geom_point(color="black", alpha=0.05, size=1) + 
  geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y, 
                                        color="red"), 
              size=2)

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

strucchange analysis

spline.data <- as.data.frame(x= spline_fit$x) |> 
  mutate(y = spline_fit$y)|> 
  rename(x = `spline_fit$x`)
plot(spline.data)

fs.data <- Fstats(spline.data$y~1) 
bp.data <- breakpoints(spline.data$y~1)
summary(bp.data)

## 
##   Optimal (m+1)-segment partition: 
## 
## Call:
## breakpoints.formula(formula = spline.data$y ~ 1)
## 
## Breakpoints at observation number:
##                       
## m = 1               78
## m = 2            70 85
## m = 3         55 70 85
## m = 4   15    55 70 85
## m = 5   15 37 55 70 85
## 
## Corresponding to breakdates:
##                                
## m = 1                      0.78
## m = 2                  0.7 0.85
## m = 3             0.55 0.7 0.85
## m = 4   0.15      0.55 0.7 0.85
## m = 5   0.15 0.37 0.55 0.7 0.85
## 
## Fit:
##                                                    
## m   0       1       2       3       4       5      
## RSS 2321698  379911  154479  140940  140332  138992
## BIC    1298    1126    1046    1046    1054    1063

plot(bp.data)

breakpoint.values=c(55,70,85); spline.data[breakpoint.values, c("x")]

## [1] 33.94309 35.04506 36.14703

segement analysis

Use out from strucchange summary to number of breakpoints as well as approximate location of breakpoints

my.lm <- lm(data=spline.data, y ~ x) 

spline.data[breakpoint.values, c("x")]

## [1] 33.94309 35.04506 36.14703

my.seg <- segmented(my.lm, 
                    seg.Z = ~x, 
                    psi=list(x=spline.data[breakpoint.values, c("x")]))
#summary(my.seg) 

my.seg$psi ; slope(my.seg)

##         Initial     Est.     St.Err
## psi1.x 33.94309 34.13619 0.04916754
## psi2.x 35.04506 34.92604 0.06043609
## psi3.x 36.14703 35.75675 0.07509539

## $x
##             Est. St.Err.  t value CI(95%).l CI(95%).u
## slope1   0.24934 0.64634  0.38578   -1.0343     1.533
## slope2  78.34700 7.65450 10.23500   63.1450    93.550
## slope3 161.85000 7.65450 21.14400  146.6400   177.050
## slope4 221.18000 2.89310 76.44900  215.4300   226.920

my.fitted <- fitted(my.seg)
my.model <- data.frame(temperature = spline.data$x, rate.output= my.fitted)

ggplot(my.model, aes(x=temperature, y=rate.output)) + geom_line()

ggplot(posterior_draws, aes(x=Temperature, y=.linpred, color=.linpred)) +
  geom_point(alpha=0.05, size=1) + 
  geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y), color="red", 
              size=3) + ylab("rate.output")+
  geom_vline(xintercept = my.seg$psi[1,2]) + 
   geom_vline(xintercept = my.seg$psi[2,2])+ 
  geom_vline(xintercept = my.seg$psi[3,2])+ 
   geom_vline(xintercept = my.seg$psi[1,1], color="blue") + 
  geom_vline(xintercept = my.seg$psi[2,1], color="blue") + 
  geom_vline(xintercept = my.seg$psi[3,1], color="blue") + 
  geom_line(data=my.model, aes(x=temperature, y=rate.output), color="black", size=2) +
  theme_classic()

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

High-latitude - 31C

posterior_draws <- lat_resp_dat2 |> 
  group_by(region, dev.temp,dev.temp.region) |> 
  data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature), 
                          to=max(lat_resp_dat2$Temperature), 
                          length.out=100), 
            mass=0.00148299) |> 
  add_linpred_draws(model1.gamk4, 
                   n=300, 
                   re_formula=NA) |>  
  unite("exp_group", c("region","dev.temp"), remove=FALSE) |> 
  filter(dev.temp.region == "31_leading")  

spline_fit <- smooth.spline(posterior_draws$Temperature, posterior_draws$.linpred, nknots = 4) 

ggplot(posterior_draws, aes(x=Temperature, y=.linpred)) + 
  geom_point(color="black", alpha=0.05, size=1) + 
  geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y, 
                                        color="red"), 
              size=2)

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

strucchange analysis

spline.data <- as.data.frame(x= spline_fit$x) |> 
  mutate(y = spline_fit$y)|> 
  rename(x = `spline_fit$x`)
plot(spline.data)

fs.data <- Fstats(spline.data$y~1) 
bp.data <- breakpoints(spline.data$y~1)
summary(bp.data)

## 
##   Optimal (m+1)-segment partition: 
## 
## Call:
## breakpoints.formula(formula = spline.data$y ~ 1)
## 
## Breakpoints at observation number:
##                       
## m = 1               78
## m = 2            70 85
## m = 3         55 70 85
## m = 4      26 55 70 85
## m = 5   17 32 55 70 85
## 
## Corresponding to breakdates:
##                                
## m = 1                      0.78
## m = 2                  0.7 0.85
## m = 3             0.55 0.7 0.85
## m = 4        0.26 0.55 0.7 0.85
## m = 5   0.17 0.32 0.55 0.7 0.85
## 
## Fit:
##                                                          
## m   0        1        2        3        4        5       
## RSS 995002.7 164741.2  67258.5  60638.9  57568.3  57272.2
## BIC   1213.5   1042.9    962.5    961.4    965.4    974.1

plot(bp.data)

breakpoint.values=c(55,70,85); spline.data[breakpoint.values, c("x")]

## [1] 33.94309 35.04506 36.14703

segement analysis

Use out from strucchange summary to number of breakpoints as well as approximate location of breakpoints

my.lm <- lm(data=spline.data, y ~ x) 

spline.data[breakpoint.values, c("x")]

## [1] 33.94309 35.04506 36.14703

my.seg <- segmented(my.lm, 
                    seg.Z = ~x, 
                    psi=list(x=spline.data[breakpoint.values, c("x")]))
#summary(my.seg) 

my.seg$psi ; slope(my.seg)

##         Initial     Est.     St.Err
## psi1.x 33.94309 33.54958 0.02768565
## psi2.x 35.04506 34.58443 0.02256263
## psi3.x 36.14703 35.58101 0.02554307

## $x
##            Est. St.Err. t value CI(95%).l CI(95%).u
## slope1  -7.8319 0.24673 -31.743   -8.3219   -7.3418
## slope2  31.0790 1.61930  19.193   27.8630   34.2960
## slope3  91.3370 1.61930  56.403   88.1200   94.5530
## slope4 138.6300 0.76778 180.560  137.1100  140.1600

my.fitted <- fitted(my.seg)
my.model <- data.frame(temperature = spline.data$x, rate.output= my.fitted)

ggplot(my.model, aes(x=temperature, y=rate.output)) + geom_line()

ggplot(posterior_draws, aes(x=Temperature, y=.linpred, color=.linpred)) +
  geom_point(alpha=0.05, size=1) + 
  geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y), color="red", 
              size=3) + ylab("rate.output")+
  geom_vline(xintercept = my.seg$psi[1,2]) + 
   geom_vline(xintercept = my.seg$psi[2,2])+ 
  geom_vline(xintercept = my.seg$psi[3,2])+ 
   geom_vline(xintercept = my.seg$psi[1,1], color="blue") + 
  geom_vline(xintercept = my.seg$psi[2,1], color="blue") + 
  geom_vline(xintercept = my.seg$psi[3,1], color="blue") + 
  geom_line(data=my.model, aes(x=temperature, y=rate.output), color="black", size=2) +
  theme_classic()

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

03. CTmax_analysis_resp_analysis

Elliott Schmidt

2024-05-27

Scenario

Load packages

Set working directory

Import data

Data manipulation

Exploratory data analysis

2nd order polynomial

3rd order polynomial

Fit model [random factors]

Priors

Models

Null

Model1

Model2

Model3

Model4

LOO

Fit model [fixed factors - polynomials]

Prior investigation

Priors

Models

Model1.p1

Model1.p1

DHARMa

MCMC diagnostics

model1.p2

model1.p2

DHARMa

MCMC diagnostics

model1.p3

model1.p3

DHARMa

MCMC diagnostics

LOO

Fit model [fixed factors - GAM]

GAM model k=4

GAM model k=4

DHARMa

MCMC diagnostics

GAM model k=3

GAM model k=3

DHARMa

MCMC diagnostics

Partial effects plots

conditional_effects

Model investigation

summary

R2

tidyMCMC

gather draws

bayesplot

emmeans - pariwise

emtrenns - pariwise

probabilities

summary figure

Piecewise linear regression

Low-latitude - 28.5C

strucchange analysis

segement analysis

Low-latitude - 30C

strucchange analysis

segement analysis

Low-latitude - 31C

strucchange analysis

segement analysis

High-latitude - 28C

strucchange analysis

segement analysis

High-latitude - 30C

strucchange analysis

segement analysis

High-latitude - 31C

strucchange analysis

segement analysis